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Abstract 

Homogeneous shear flows (with constant strainrate dvx/dy) are generated with the Doll's and 
Sllod algorithms and compared to corresponding inhomogeneous boundary-driven flows. We use 
one-, two-, and three-dimensional smooth-particle weight functions for computing instantaneous 
spatial averages. The nonlinear normal stress differences are small, but significant, in both two 
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■ and three space dimensions. In homogeneous systems the sign and magnitude of the shear- 

plane stress difference, Pxx — Pyy, depend on both the thermostat type and the chosen shearflow 
algorithm. The Doll's and Sllod algorithms predict opposite signs for this normal stress differ- 
ence, with the Sllod approach deflnitely wrong, but somewhat closer to the (boundary-driven) 
truth. Neither of the homogeneous shear algorithms predicts the correct ordering of the kinetic 
temperatures: Txx > Tzz > Tyy. 

PACS numbers: 02.70.Ns, 45.10.-b, 46.15.-x, 47.ll.Mn, 83.10.Ff 
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I. INTRODUCTION 

In the present work, we use nonequilibrium molecular dynamics [l| to study microscopic 
simulations of "simple shear flow" (also called "plane Couette flow"): 

v^oc y-^P^y = Py^ = -r][{dvcc/dy) + (dvy/dx)] = -r]{dv^/ dy) = -r]e . 

The flow is in the x direction so that the tensor Vf has only one nonzero element, (Vf = 
dvx/dy = e. We use the symbol P for the (symmetric second-rank Cauchy) pressure tensor 
(positive in compression and the negative of the stress tensor, P = —a); v for the (time- 
and space-dependent) hydrodynamic flow velocity; rj for the shear viscosity; and e for the 
magnitude of the imposed (or measured) shear strain rate. The various simulation types 
we consider here were designed to clarify the relationships between periodic homogeneous 
methods, thermostated everywhere, and flows with moving thermostated boundaries. All 
the methods we use are consistent with Green and Kubo's linear response theory at small 
rates of shear 21] • We are specially interested in characterizing and understanding the 
nonlinear shearplane stress difference: [ Pxx — Pyy = o'yy — Oxx ] which arises in sufficiently 
small systems at sufficiently large shear rates. In carrying out microscopic simulations 
both the boundary conditions and the thermostats or ergostats which control the flow 



need to be carefully considered 

Though simple shear flow is "stationary", fluctuations in local properties necessitate 
averaging, both in time and in space. Here we reduce the importance of these fluctuations 
by usiug spatial averag,.g techniques bonowed from smooth-part.ele eo.t.uuum Simula- 

tion methods[3]. We measure instantaneous spatially- averaged How velocity, temperature, 
and pressure-tensor components. 



The two best-known homogeneous microscopic methods, the Doll's tensor[4i| and Sllod 
algorithms [5], treat fluid rotation differently, leading to qualitatively different predictions 
for the nonlinear normal stress difference. Boundary-driven flows can help to resolve this 
disagreement Useful flows need to satisfy four conditions: the spatial scale L 

of these flows needs to be large enough (relative to the particle size), but not too large 
(to avoid turbulence), with flow velocities v large enough (to emerge above fluctuations), 
but not too large (again, to avoid turbulence), in order to provide useful information. 
The relatively greater importance of fluctuations in two dimensions is responsible for 
the reduced utility of viscosity there, as shown in Fig. 1. Additionally, the frictional 
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boundary-fluid interaction needs to be sufficiently strong to prevent excessive boundary 
slip. 

In the present work, we carry out both two- and three-dimensional simulations using all 
three approaches (Doll's, Shod, and boundary-driven) for two simple, and rather similar, 
pairwise-additive repulsive potentials. We focus here on the normal stress difference, 
Pxx ~ Pyy = o'yy — Oxx^ in the shearplane. We will see that the Doll's and Shod algorithms 
typically predict different signs for the simple monatomic dense fluids considered here. 
Both these homogeneous algorithms yield clearcut results, insensitive to system size. 

In the ea..,v da., of such .hea.flow ..u.atio.Q it was thought that the two- 
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dimensional viscosity might depend logarithmically on system size[ll|]. More recently, 

Green-Kubo stress-stress correlation 



careful studies 
function: 

/■oo 

ri = {V/kT) / ( P^y{Q)P.,y{t) ),^dt , 

Jo 

indicate no such dependence in dense fluids. The work we carry out here is consistent 
with this lack of size dependence. Here we find that the coefficient of the hypothetical 
logarithmic viscosity contribution can be no larger than 10~^ in the natural reduced units 
of atomic size, mass, and velocity. The more realistic boundary- driven flows necessarily 
entail larger fluctuations and considerable size dependence. We are nevertheless able to 
determine the sign and the size of the normal stress difference for such flows so as to 
characterize the errors inherent in the homogeneous algorithms. 

This paper is organized as follows: in Sec. II we lay out the macroscopic description 
of the problem; in Sec. Ill the microscopic description; In Sec. IV we describe the 
two homogeneous algorithms emphasizing their similarities and differences; in Sec. V 
we describe the boundary-driven algorithm used in the present work; Sec. VI describes 
the smooth-particle spatial averaging method used in analyzing results from simulations; 
Sec. VII describes numerical results for two similar (smooth repulsive) forcelaw models 
in two very- different density regimes and in two space dimensions; Sec. VIII describes 
corresponding three-dimensional results for one of these forcelaw models; Sec. IX lists the 
conclusions we have drawn from this work. 
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Figure 1 
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Figure 1: Limits imposed on system size L and boundary velocity v by atomistic size, thermal 
fluctuations, turbulence, and Shockwaves. Viscosity is a useful concept for flows in the enclosed 
areas. The figure here is constructed for a fluid at unit mass, number density, temperature, 
viscosity, and heat conductivity. L is measured in units of the microscopic particle size and v is 
measured in units of the thermal velocity. Modeled after Ref . [6] . 

II. MACROSCOPIC DESCRIPTION OF SIMPLE SHEAR FLOW 



Classical fluid flow can be modeled and understood from either the microscopic or 
the macroscopic standpoint. In the microscopic description individual particles obey the 
ordinary differential equations of motion of nonequilibrium molecular dynamics 

{mf = mv = FA + FB + Fc + FD}, 

and everything follows from the functional forms of the assumed atomistic, boundary, 
constraint, and driving forcesQ, In molecular dynamics, just as in continuum hy- 
drodyamics, the pressure tensor is (defined to be) the comoving momentum flux, comoving 
relative to the flow velocity. In a many-body system with forces derived from the pairwise- 
additive pair potential 0(r), the pressure tensor is made up of {i,j} pair contributions as 
well as individual particle convective contributions 



i<j i 



X'lj — Xi Xj , t/ij — t/i t/j , Tjj — |rj I , Fij — ^ i^j^iXij) ■ 

It is evident that the pair-potential pressure tensor is symmetric, with P^y and Py^ equal. 

In the macroscopic description the continuum field variables [ such as the mass density 
p{r,t), the velocity v{r,t), and the energy per unit mass e{r,t) ] obey partial differential 



equations: 

p = -pV • V ; 

pr — pv — —V • P = V • (J ; 

pe = -Vv : P -V -Q . 

In these general continuum field equations the constitutive relations for the stress tensor 
a = —P and the heat flux vector Q distinguish one material from another. In both 
approaches a computational algorithm for solving the equations is needed. Its imple- 
mentation gives { r(t),v{t) } in the microscopic case and { p{r,t),v{r,t),e{r,t) } in the 
macroscopic case) . In either case, a well-posed problem also requires boundary and initial 
conditions, constraints, and driving forces. 

First consider the simplest model system illustrating stationary simple shear: imagine 
an incompressible Newtonian fluid, with constant shear viscosity rj, and which also follows 
Fourier's linear heat transport law with a constant heat conductivity k.: 

Pxy = -r][ {dvx/dy) + [dvy/dx] ] ; 
pe = kW^T -P:Wv . 

We ignore thermal expansion, so that the mass density p is constant. We denote the local 
thermodynamic variables in the conventional way: temperature T{r), pressure tensor 
P(r), and internal energy per unit mass e(r). We adopt the colon convention in the 
tensor product A : B to indicate a sum over all four A^jB^j terms in two dimensions, and 
all nine such terms in three dimensions. In some texts the alternative sum (immaterial 
for the symmetric tensors considered here) AijBji is used. 

Simple shear flow for this bare bones textbook model is perhaps the simplest imaginable 
nonequilibrium flow problem. It gives a linear variation of velocity in space along with a 
quadratic variation of temperature. Simple shear flow can be driven by two moving parallel 
boundaries, both of them at temperature Tb, and able to absorb heat and to impose their 
boundary velocities, Vx = ^v, on a two-dimensional strip or three-dimensional slab of 
model fluid of thickness L: 

Vx{y = ±L/2) = ey . 

The stationary macroscopic description of such a flow, for the model Newtonian fluid 
with Fourier heat conduction, has a constant stress tensor, a linear velocity proflle and 
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Figure 2 
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Figure 2: The simple shear flow indicated in the center panel corresponds to the linear velocity 
profile (shown to the left) and quadratic temperature profile (shown at the right). The nonlinear 
effects considered in this paper give relatively small deviations from these idealized profiles. 

(because V^T is constant) a quadratic temperature profile: 

(dvx/dy) = i = 2v/L = -Pxy/v ; 
AT{y) = T{y) - = {r^e^ /2k)[{L/2) - y][{L/2) + y] . 
The maximum temperature difference, relative to the boundaries' temperature Tb, 
AT{y) < AT^,, = Ar(0) = r^e^LVS/t = r]vy2K , 

occurs at the mi dp lane y = 0. See Fig. 2 for a schematic illustration of this prototypical 
simple shear flow. 

This stationary solution satisfies energy balance, with the rate at which heat is gener- 
ated throughout the volume V (necessarily the same as the rate at which external work 
is done, —W) equal to the rate at which heat is transferred through the two boundary 
walls of area (length in two dimensions) A, at y = ±L/2: 

-W = -P^yVi = 7]Ve^ = KA[{dT/dy).L/2 - {dT/dy)+L/2] • 

This macroscopic description of shear flow guides our interpretation of microscopic simu- 
lations. We focus on the complications caused by fluctuations and nonlinearities in what 
follows. 



III. MICROSCOPIC DESCRIPTION OF SIMPLE SHEAR FLOW 



Molecular dynamics was developed over 50 years ago 17|, ll8|, ll9| and soon gave rise to 
successful interpretations of equilibrium properties based on hard-sphere perturbation- 
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theory 



20] analogous to Enskog's hard-sphere understanding of dense-fluid transport 



properties 211] . Fie 



conducting flows 



24 



-driven diffusion 22!], shear and bulk viscous flows 



231], and heat 



251 ] all came to be simulated with a variety of algorithms. Spe- 



cial boundary conditions and computational thermostats were developed to model these 
flows. Throughout this development the limiting case of Green and Kubo's linear theory 



of transport served ClS db guide 



26|. 



In the past thirty years simulation has come a long way. Billions of atoms can be 



simulated now 27|. Nonequilibrium simple-shear algorithms have been developed for di- 



atomic and polyatomic molecules 



28 



291 ]. not just simple fluids. More recently, alterna- 



tive shear-flow algorithms have been developed for periodic irrotational "elongational" 



flows 



31 



321 ]. flows with a steady stretching in one direction and a simultaneous 



shrinking in a perpendicular direction. It might be thought (as it once was 33|]) that such 
a flow could not be followed forever, but a clever choice of periodic boundaries makes 
it possible to study steady thermostated elongational flows. The conce ptu al difficulties 
involved in deriving these algorithms have led to a spirited literature 3^, |35|, l36| as to the 



"correctness" of the various algorithms. Such discussions can easily lead outsiders to the 



impression that it is hard to distinguish "correct" from "incorrect" a 



35 



'he controversial aspects of these shear-flow algorithms 



28 



29 



gorithms 



30 



31 



32 



33 



34 



361 ] led us to reconsider the problem. It seemed to us that a fresh look at the basic 



algorithms for monatomic simple-shear flows would help to develop a perspective clarifying 
this situation. It is evident that the nonlinear aspects of the computer algorithms are to 
some extent arbitrary, as the only true guidelines for correctness are consistency with the 
well-known and well-accepted linear flow theory described by Newtonian viscosity and 
Fourier heat conduction. 

In viscous shear flow, which we consider here, any reasonable algorithm needs to satisfy 
the requirement that the shear viscosity for small strain rates agrees with Green and 
Kubo's linear-response relationj^, 26] linking the viscosity to equilibrium fluctuations in 



the shear stress a^y = —Pxy- 



ri = {V/kT) / (P.,(0)P,,(t))eqrft 



- p 

-f xy 



' e ' dy ^ dx 

For simplicity we consider systems with pairwise-additive forces {Fij} derived from a 
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potential function $ = ^ 0: 



{F,, = -V.0(|r.-r,|) = -F,, }, 
In such a system both the energy E, 

E = $ + = ^ 0,,. + ^ p^/2m , 

i<j i 

and the microscopic pressure tensor P, 

i<j i 

are sums of two-particle potential and single-particle kinetic contributions. 

There are two interesting ambiguities in the definition of pressure in a particulate 
system. The contributions of the potential pair interaction terms { {Fr)ij } to PV 
need to be allocated spatially. Either delta-function contributions at and rj, or at 
{vi + rj)/2, or smoothed distributions centered at these locations can be used. The old 
"Irving-Kirkwood" preference for delta functions is motivated more by analytic conve- 
nience than by any physical considerations. A smoothed approach is certainly preferable 
in computational work. 

The kinetic part of PV is ambiguous too. How is the "comoving" velocity to be 
defined in a system with transient velocity fluctuations? Either an instantaneous or a 
time-averaged hydrodynamic velocity can be chosen. Only the kinetic part of PxyV is 
affected by this choice. The simplest versions of the two choices are as follows: 

Pxyi^)^ = X^^^bx ~ '^xit)][Vy — Vy{t)] (instantaneous) 

i 

or 

Pxy{i)y = - {vx)time\[vl " {vy)timc] (time avcragcd) . 

i 

The kinetic part of the pressure tensor is the usual definition of the instantaneous ki- 



netic temperature for N particles 37|]. In three space dimensions the usual definitions of 
\TxxiTyy^Tzz\ are. 



N 

2 



N 

NkTyyit) = Y,HVI- {Vy{t))f : 
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N 

NkT^,{t) = Y,m{vi-{vM)f , 

where the average velocity at time t is computed by summing individual particle 
velocities. 

The instantaneous hydrodynamic velocity {vx{t),Vy{t)) has to be estimated numeri- 
cally. The simplest way to do this in a system composed of moving particles, is to use a 
weighted (smoothed-particle) average of nearby particle velocities: 



!(t) ='^ViWir/'^Wir ; Wir = W{\r - Til) 



where the time-independent scalar weight function Wir has a finite spatial range, and 
emphasizes the contributions of those nearby particles to the location r where v{t) is to be 
computed. The details of this spatial average are described in Section VII. Fortunately, the 
main conclusions of this work are independent of the necessarily arbitrary choice between 
using instantaneous and time-averaged flow velocities. By choosing to study stationary 
simple shear we avoid the additional complexities associated with time averaging flows 
which have an explicit time dependence in their boundary conditions. 

The atomistic simulation approach is also complicated by heat and by fluctuations. 
Because viscosity is dissipative, a viscous shear flow necessarily generates heat. Particle 
motions and collisions in systems tractable with molecular dynamics necessarily exhibit 



fluctuations in all their local properties 38j. Despite the time-reversibility of the un- 



derlying equations of motion, the overall longtime-averaged character of these flows is 
necessarily dissipative, in accord with macroscopic hydrodynamics and the Second Law 



of Thermodynamics 38|, l39|, |40|, l41 |. 



To simulate a steady viscous flow, the heat generated needs to be removed. This can be 
done with any one of many schemes, all of which are based on time-reversible constraint 
forcesQ, 3]. Such constraint forces can play the role of a feedback-based thermostat or 
ergostat, 

{ -^Constraint = -(p }^k = OoTE = 0. 

Typical constraint-force choices keep the kinetic energy K or the total energy E flxed, or 
allow these energies to fluctuate about a specifled mean value in a way consistent with 



Gibbs' statistical mechanics at equilibrium 42, |43|]. For simplicity, we restrict ourselves 



here to "Gaussian thermostats" [ so named after Gauss' Principle of Least Constraint 



44l | which flx the kinetic energy K or the total energy ii^ of a particular set of degrees 

9 



of freedom by imposing feedback-based constraint forces, { -^constraint = —CP }■ The 
Gaussian "friction coefficient" ( constrains the momenta { p } contributing to K or E. 
We next describe the two best-known computational algorithms for simulating simple 
shear. 



IV. DOLL'S AND SLLOD ALGORITHMS FOR SIMPLE SHEAR FLOW 



Shear flows were the first application of homogeneous nonequilibrium molecular 
dynamics^]. Steady inhomogeneous Shockwaves 45] and boundary-driven shear fiows0, 
7, y, had both been simulated earlier. Special boundary conditions, or external forces, 
as well as thermostats had to be developed for the homogeneous shear fiows. The nu- 
merical applications of these ideas, in the early 1970s, preceded their formal theoretical 
development by decades. The Doll's Tensor Hamiltonian for viscous fiows was discovered 
in 1985 4| ; Dettmann and Morriss' Hamiltonians (for the Nose-Hoover and Gauss ther- 
mostats) were discovered in 1996 [4^ 4^]; and a proper formulation of elongational fiows 
first appeared in 1998 jsot. 

Simple homogeneous shear fiow, with the x velocity proportional to the y coordinate, 

e = dvxidy *~ 



Vv 



dvx 


dvy_ 










dx 


dx 




dvx 


dVy_ 




e 





dy 


dy 







48 



developed 



can be implemented with periodic "Lees-Edwards" boundary conditions 
independently by Bill Ashurst in his Ph D thesis work ] . See page 26 of Ref . [6] for a 
brief description of Ashurst 's algorithm. The corresponding fiow is illustrated in Figure 
3. 

This simple shear fiow, with ^ nonzero, is "rotational" , in the sense that the clockwise 
rotation rate, — u;, is nonzero, and equal to half the strain rate: 

-u = [{dv^/dy) - {dvy/dx)]/2 = e/2 . 

In our specific special case Vf can be written as a sum of irrotational and rotational 
contributions: 
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Figure 3 
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Figure 3: Periodic 100-particle shear flow, showing eight periodic L x L images of the central 

10 X 10 box. The particle velocities in the images above/below the central A^-particle box differ 
from those of the central box by iteL. 



See Figure 3 for an illustration of the periodic boundary conditions consistent with such 
a flow. A stationary nonequilibrium state can result if the irreversible heating rate, 
—W = rje^V is appropriately compensated by momentum-dependent thermostat forces. 

The Doll's and Sllod algorithms are straightforward possibilities for simulating such a 
flow, and the differences between them are small (of order e^) at moderate strain rates. 
Ashurst's periodic shear algorithm preceded the formal derivations of differential evolution 
equations for the shear flow and its associated thermostats. His flnite-difference algorithm 
was equivalent, for small timesteps, to a thermostated version of what is now called the 
Sllod algorithm: 



{ i = (Px/m) + ey ; y= (Py/m) ; 



Fx - ePy - CPx 



Py 



CPy } ■ 



In the absence of the thermostat forces { —(p }, these motion equations give an energy 
change exactly consistent with the rate at which thermodynamic work is performed by 



the instantaneous shear stress —P. 



xy 



In the alternative Doll's- Tensor approachjj], a coordinate- momentum {q,p) sum [ giv- 
ing rise to the "Kewpie-DoU" name ] is added to the usual Hamiltonian. The added term 
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includes the strain rate e: 



"^^Doll = "^^Usual + e • 

i 

The resulting Hamiltonian equations of motion, 

dp ' dq ' 

are only slightly different to the closely related "Sllod" equations. The Doll's tensor 
motion equations are 

{x = (px/m) + ey ; y = (py/m) ; p^ = ; py = Fy - ep^ } . 

In the absence of heat-extracting thermostats, both approaches, Sllod and Doll's, provide 
an energy change exactly consistent with thermodynamics: 

HDoii(g,P, Vt;) = ^ 7iusuai(g,p) = -{d/dt)^qp : Vf = -ePaiyV^ • 

The extra momentum-dependent forces, {—ipy} for Sllod and {— epx} for Doll's, model 
rotation. In a rotational flow in the xy plane the rotation rate a; is given by 



1 

^"2 



dvy dvx 



dx dy 

and corotating momenta would include rotational contributions: 

Px = -^Py ; Py = • 

The Coriolis accelerations corresponding to rotation are similar in form to the momentum- 
dependent forces included in the Doll's and Sllod algorithms. 

It should be noted that the ambiguity in treating the atomistic rotation associated with 



shear has a relatively well-known analog in continuum mechanics [49|]. In the plastic flow of 
solids the total macroscopic strain, idealized as a sum of elastic and plastic contributions, 
can only be computed by time integration: J edt e. Because stress depends on strain, 
calculations involving plastic flow require the simultaneous time integration of an equation 
for stress evolution, J &dt —>■ a. If the effect of rotation on the stress is included (imagining 
that the stress rotates as if it were embedded in a rigid body) the "Jaumann stresses" 
result. These laboratory- frame stresses satisfy the evolution equations: 

^xx 2tUCT2,y , ^yy ~\~'^^(^xy 1 ^xy ^Ip'xx ^yy^ ■ 
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A more complicated boundary-driven shear flow can be driven by imposing appropriate 
boundary conditions, including external thermostats or ergostats, on a Newtonian region 
with the usual equations of motion: 



{x= (p^/m) ; y = (py/m) ; p^ = ; py = Fy } . 

The details of our boundary- driven simulations are given in Sec. VI. In his thesis work 
in the early 1970s Ashurst 7|] found consistent shear viscosities for both boundary-driven 
and homogeneous periodic shear flows. More recently Liem, Brown, and Clarke made an 
effort to characterize the very nonlinear effects we seek in the present work in shear flow. 
They compared homogeneous shear (presumably using the Sllod approach) with a very 

n 

large scale boundary-driven flow [9]. Unfortunately the limitations inherent in using such 
large systems made the nonlinear effects too small to measure. 



V. BOUNDARY-DRIVEN SHEAR FLOW 



To avoid dealing with free surfaces we have chosen to simulate fully-periodic four- 
chamber systems of the type shown in Figure 4. Posch and Hoover used the same geometry 
in an investigation of heat flowjsol. In two dimensions the system is made up of four Lx L 
cells, each with an area one-fourth the total, = N/4. To eliminate surface effects the 
system is periodic in both x and y, with two steadily-moving chambers thermostated at 
a fixed kinetic temperature. The motion is induced by using steadily moving tethers. 
The remaining particles obey Newton's motion equations. This four-chamber geometry 
has a statistical advantage. The two Newtonian regions provide independent estimates of 
the stress and temperature tensor profiles. The differences between the two Newtonian 
regions provides a useful indication of the estimates' reliability. See, for instance, the 
normal-stress profiles for the two regions shown in Figs. 8, 10, and 11. 

Though we do not give the results here, we have also considered a two-chamber periodi- 
cally shearing flow geometry, in both two and three space dimensions. In the two-chamber 
approach the moving tethers' x velocities are 

{ ±UL , ±3eL ; ±5eL ; . . .} , 

in the odd-numbered cells, where L is the cell width, while the motion in the even- 
numbered cells is purely Newtonian. Just as in the four-chamber case periodic shearing 
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Figure 4 
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Figure 4: Four-chamber L x 4L periodic flow with 4 x 400 particles. Chamber 1 (filled circles 
at the bottom) moves to the right at speed +eL/2; chamber 3 (filled circles, third chamber from 
the bottom) moves to the left at speed —eL/2. This imposes the nominal strain rates on the 
two Newtonian chambers 2 and 4 (open circles). Periodic boundaries apply in both the x and 
the y directions. Note the ordering effect of the moving tethers in Chambers 1 and 3. 

boundaries are used in all the space dimensions. The results from the two-chamber and 
four-chamber flows are not significantly different and are therefore not elaborated here. 

The multi-chamber shear flows are induced by tethering the boundary particles to 
steadily moving lattice sites (a moving square lattice or simple cubic lattice, for con- 
venience) with a simple quartic tethering potential, just as in the 0*^ model for heat 



conduction 



51 



52|: 



4[r-ro(t)r ; ro(t) = (±eL/2, 0) 



Here the sites move. 100 is a good choice for the force constant ^4 and gives better 
equilibration than tethers with alternative power-law exponents. 
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In order to prevent diffusive mixing between the Newtonian particles and the bound- 
ary particles it is useful to include a smooth repulsive boundary force returning errant 
Newtonian particles toward their chambers. For convenience this force has the same form 
as the boundary tethering forces used in the moving chambers. The fixed speed of the 
tethers, |ro(^)|; is chosen to impose overall strain rates ±e on the two Newtonian regions. 
Whether or not the strain rate actually penetrates into the Newtonian regions depends 
upon whether or not "slip" occurs at the boundary, as we discuss in Sec. VIIIB. 

Because the heat generated by the shear varies as in D dimensions, but must flow 
across a boundary of area (or length) L^~^ something "interesting" necessarily occurs 
as L increases. Eventually the driving chambers must become decoupled from the bulk 
Newtonian chambers so that both the work done and the corresponding heat generated 
actually decrease with increasing strain rate. In favorable cases, with good frictional 
coupling at the boundaries, this driving technique provides good velocity profiles and 
allows nonlinear effects to be determined. We turn next to the instantaneous smooth- 
particle spatial-averaging procedure required for analyzing the simulations. 



VI. SMOOTH-PARTICLE SPATIAL AVERAGES 

The spatial averaging algorithms our simulations require are borrowed from a contin- 
uum technique, "SPAM" jsl. Smooth particle applied mechanics ("SPAM") is a technique 
for solving the partial differential continuum equations (continuity, motion, and energy) 
for the evolution of the density, velocity, and energy. It makes use of a normalized weight 
function with a maximum range /i, w{\r\ < h). The weight function describes the spatial 
extent of a representative particle of mass m. The weight function is formulated with at 
least two continuous derivatives everywhere in order that the first and second continuum 
derivatives of smooth particle sums [ corresponding to instantaneous local quantities such 
as Vp(r, t), Vf(r, t), and V^T(r, t) ], are everywhere continuous in both space and time. 
This continuity and differentiability facilitates a smooth transition between particle and 
continuum analyses. 

Consider the two simplest cases, mass and momentum sums. The density p{r), and the 
momentum density p{r)v{r) at the location r are definedhy summing up the contributions 
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of all particles {j} within the maximum range h of that location: 

p(r) = m w{\r — Tjl) = m ■ 



J 



p{r)v{r) = m Vjw{\r — |) = m WrjVj . 
j 

An advantage of this formulation is that the continuum continuity equation, 

is satisfied exactly by the interpolated smooth particle fields: 

^P('^)\ sr^ ( ^'^^'^ ~ 3\) 

I . I ^ 

J ^ ^ i 

Notice that the spatial gradient operator, Vr, affects only w{r — rj) and not the individual 
point-particle properties {vj}. 

By choosing properly normalized weight functions spatial averages can be computed 
in one, two, or three spatial dimensions. The pressure tensor in the vicinity of a particle 
or at a grid point in three dimensions would be computed by using the t/iree-dimensional 
weight function: 

wsoir) = (105/167r/i^)[l - 6{r/hy + S{r/hf - 3(r//i)^] / A-nr'^W3D{,r)dr = 1 . 



This particular weight function was discovered and used by Lucy in 1977[53|. It is the 
simplest normalized polynomial with (1) a maximum at r = 0; (2) a finite range r < h, 
and (3) two continuous derivatives vanishing a.t r = h. The corresponding one- and two- 
dimensional weight functions (for averages in thin strips or slabs, and for averages at 
particles or grid points in two-dimensional problems, are 

W2D{r) = {5/nh'^)[l-Q{r/hf + 8{r/hf-3{r/h)'^]^ / 27Trw2D{r)dr = 1 . 

Jo 

eh 

\2 I o/^/^,^3 Q/^/L\4l 



Jo 

The smooth-particle equations of motion for the time development of the individual par- 
ticle velocities {vj} are generally formulated so as to conserve linear momentum exactly. 
The failure of this approach to conserve angular momentum is a relatively subtle point 



worthy of more research investigation 54 ] . 
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VII. NUMERICAL RESULTS: TWO DIMENSIONS 



A. Homogeneous Simulation Results with Lucy's Potential 

work, reduce numerical integration errors, and to make 



To simplify the numerica 
contact with existing data 55 



56 



57j, we initially chose to use Lucy's pair potential in 
two space dimensions. This pairwise-additive potential is identical to the weight function 
'W2D{r) of the previous Sec. VI: 

^ ■[! - + Sx^ - 3x^] ; X = ^ < 1. 



' 7r/i2^ " h 

We can make a rough estimate of the potential contribution to the energy and pressure 
for this potential by assuming a random distribution of particles. Viewed as a statistical- 
mechanical potential function for molecular dynamics in two dimensions, Lucy's potential 
(due to the normalization of the Lucy weight function) then corresponds to the simple 
equation of state for a two-dimensional ideal gas: 

^ 2V 2V 

The energy expression follows from the normalization condition, while the pressure ex- 
pression follows from the virial theorem: 

t2 rh J^2 



1 T\T2 rn 



2V 

At unit mass and number density and at a strain rate of e = 0.05 the shear viscosity and 



normal stresses (for Sllod) were measured precisely in 1995 55|, |56|]. Because the results 
are insensitive to the number of particles used we list in Table I below only a single 
representative set of simulations for = 400. The insensitivity of the viscosity to the 



algorithm type is also evident in Daivis recent work [58 1. 

We chose h = 3, for which the Sllod algorithm shear viscosity (with a single thermostat 
variable) has previously been computed. At unit density, with L x L = N, the average 
number of pa,r mteractions approximately UN. We ako used exactly the same Lucy's 
function as a weight function[3] for carrying out spatial averages. Table I lists the pressure- 
tensor components using both the Doll's and Sllod algorithms. The first four simulations 
are "isothermal" with both the average temperature, [T^^ ~^Tyy)/2 thermostated (SI and 
Dl), and with both temperatures separately thermostated (S2 and D2), by using two 
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control variables, C,x and C,y. The last two simulations (SE and DE) constrain the internal 
energy rather than the temperature. Notice that in all these simulations the average 
pressure is fairly close to the simple virial-theorem estimate, P^V/N ~ 0.5. 
The normal stress difference is sensitive to the single-thermostat type: 

Pxx - Pyy = +0.015 (Shod) ; - Pyy = "0.017 (DoU's) , 

and is almost entirely kinetic. At the high effective density of these simulations the 
potential portion of the stress distribution is nearly isotropic. These results are essentially 
unchanged if energy, rather than average temperature, is controlled. On the other hand, 
this striking normal-stress effect disappears completely (because it is kinetic in nature) if 
the two temperatures are separately thermostated. We conclude from these simulations 
that a definitive determination of the normal stress difference requires a more realistic 
boundary-driven flow. One can have no confidence in nonlinear kinetic effects when the 
evolving kinetic energy itself can be dominated by the choice of homogeneous thermostats, 
one or two. 

Table I. Relatively high-density Lucy viscosity in two dimensions. The mass and 
number densities are unity: p = niN/V = N/V = 1. The range of the Lucy potential 
is /i = 3, so that each particle interacts simultaneously with approximately 30 neighbors. 
Space-and-time- averaged pressure tensors are given here for homogeneous Sllod and Doll's 
algorithms with a square periodic 20 x 20 = 400-particle cell with a strain rate of e = 
dvx/dy — 0.05. The kinetic temperature is fixed, kT — 0.07, using a single control variable 
C (in the two runs SI and Dl) and two separate control variables in the x and y directions, 

and Cy (in the two runs S2 and D2). The total energy is fixed at 0.500 for the runs 
SE and DE. The pressure-tensor components are given in the order xx, yy, xy with the 
kinetic, potential, and total terms indicated. The boundary conditions are periodic, with 
a total run time of 200 x 10, 000 timesteps. The fourth-order Runge-Kutta timestep is 
0.005. The average potential energies per particle for the first four runs are all in the 
range 0.4272 <^/N < 0.4274 so that the total energy in these runs is 0.497. 
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Figure 5 

0.3 
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-0.1 
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Figure 5: 400-paxticle Lucy velocity profile. See Fig. 4 for a 1600-particle analog. The 
two moving chambers correspond to the maximum and the minimum velocities. Here the two 
Newtonian chambers are so poorly coupled to the moving chambers that their mean velocities 
are nearly zero. 



Run Type 


pK pH 


pK p$ pS 
yy yy yy 


pK pf> pE 

xy xy xy 


400S1 
400D1 
400S2 
400D2 
400SE 
400DE 


0.078 0.484 0.562 
0.060 0.486 0.546 
0.070 0.485 0.555 
0.070 0.485 0.555 
0.081 0.484 0.564 
0.062 0.486 0.548 


0.062 0.485 0.547 
0.079 0.483 0.563 
0.070 0.484 0.554 
0.070 0.484 0.554 
0.064 0.485 0.549 
0.083 0.483 0.566 


-0.022 +0.002 -0.021 
-0.021 +0.001 -0.020 
-0.025 +0.002 -0.023 
-0.025 +0.002 -0.023 
-0.023 +0.002 -0.022 
-0.022 +0.001 -0.021 



B. Boundary-Driven Simulation Results with Lucy's Potential 

Because the homogeneous investigations with the Lucy potential were sensitive to 
thermostat type, we next carried out several four-chamber simulations with two boundary 
chambers moving oppositely, as shown in Fig. 4. The density and strain rate were chosen 
to match the data in Table I. These simulations produced no useful normal-stress results! 
This failure reflects the nearly negligible coupling between the two moving chambers 
and the two Newtonian chambers. See the typical velocity profile in Fig. 5. Particles 
in the two Newtonian chambers simply rest quietly between the two rapidly moving 
walls. Evidently this very dense repulsive fluid with relatively weak coUisonal forces has 
insufficient friction for boundary driving to reach strain rates with significant nonlinear 
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Figure 6 
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Figure 6: Comparison of Lucy's potential and the short-ranged soft-sphere repulsive potential 
4>. The ordinate and the abscissa are divided by their maximum values, so that both scales vary 
from to 1. 

stress differences. This same difficulty was found by Liem, Brown, and Clarke in their 
three-dimensional Lennard- Jones simulations!^. 

On the other hand, in Ashurst's thesis work his "fluid- wall" boundary driving regions 
(velocity constraints, but no tethers, were used) produced good linear velocity profiles 

. Because those simulations cor- 
respond to a much lower density (with about three interactions per Lucy particle rather 
than thirty) and much more violent collisions than those of Table I we abandoned the 
high-density Lucy simulations and took up instead soft-disk and soft-sphere simulations 
at conditions more closely resembling those of Ashurst. The forcelaw change was moti- 



vated by the desire to check our results with those from previous simulations 
new simulations are discussed in the following two subsections. 



57|. The 



C. Homogeneous Simulation Results with a Soft Disk Potential 



Here we consider the short-ranged smooth soft-disk repulsive potential 57|] , 

0(r < 1) = 100(1 - r'^Y ; r"^ = + y"^ , 

very similar in general shape to Lucy's, but with three vanishing derivatives rather than 
just two, at its maximum range of unity. The two potentials are compared in Fig. 6. 
In comparison simulations we found that the two potentials provide quite similar normal 
stresses at corresponding temperatures and densities. Here we choose unit density and 
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energy per particle. We will see that boundary- driven simulations with this potential 
choice provide good velocity profiles, as did Ashurst's similar "fluid walls" in 1972, and 
also allow comparisons with previous viscosity simulations in the same thermodynamic 
state. First we consider again Sllod and Doll's simulations with homogeneous shear. 

Just as in the Lucy simulations, the soft-disk normal stresses are quite different for 
the Doll's and Sllod algorithms and are even less sensitive to system size. Now the short- 
ranged potential contribution to shear stress is actually dominant. The results listed in 
Table II, all using a single friction coefficient, fixing the energy per particle, E/N = 1, 
show that now the potential contribution to the stress difference is the same order as the 
kinetic contribution. The Sllod pressure difference is considerably smaller in magnitude: 

Pxx - Pyy = +0.008 (Sllod) ; P,,. - Pyy = -0.054 (Doll's) , 

Because the rotated normal stress difference is equivalent to a shear stress: 

{Pxx - Pyy) 45° p 

2 ^xy : 

an alternative phenomenological description of the normal-stress effect corresponds to a 
rotation of the principal stress direction (the direction of maximum tension in the traceless 
stress tensor). The Doll's algorithm gives a clockwise shear-stress rotation while the Sllod 
rotation is counterclockwise. 

Table II. Moderate-density soft-disk viscosities. The range of the potential is unity 
so that each particle interacts, on the average, with only three neighbors at unit density. 
Space-and-time-averaged pressure tensors for homogeneous Sllod and Doll's algorithms 
using the soft-disk pair potential illustrated in Fig. 6, = 100(1 — r^)^. The energy and 
density are fixed, and equal to unity, and the strain rate, e = dv^/dy is 0.50. The pressure- 
tensor components, are again given in the order xx, yy, xy with the kinetic, potential, and 
total terms indicated. The boundary conditions are periodic, with a total run time of 
200x10,000 timesteps for N = 64, 256, and 1024; 200x2000 timesteps for N = 4096, and 
200x500 timesteps for N = 16,384. The fourth-order Runge-Kutta timestep is 0.005 and 
a single isoenergetic friction coefficient is used, as explained in Sec. II of the text. 



21 



N Type 
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The SUod program used to generate the data in Table II reproduced earlier Sllod 
results 57 1 for strain rates of 0.10 and 0.25 very well. Gass' Enskog-theory prediction {59! 
of the (linear) viscosity for the conditions of Table II, r/Enskog = 1-5, is just slightly higher 
than the value f] = 1.25 given by the data in this Table. Just as in the work of Refs. 
[10]- [12] there is no indication of any logarithmic A^-dependence in these results. Let us 
now compare these homogeneous results to those from boundary-driven simulations. 



D. Boundary-Driven Results with the Soft-Disk Potential 

With the soft-disk potential 0(r < 1) = 100(1 — r^)^ a useful boundary-driven flow 
does result if the strain rate is moderate and the system is not too large. The boundary- 
driven flows are sensitive to system size. This is because the heat generated in a shear 
flow, of order for L large eventually overwhelms the capacity of the boundary, of order 
L^~^, to absorb it. The maximum strain rate that can be reached by boundary- driven 
flows is therefore limited by the heat conductivity as well as the efficiency of the frictional 
thermal contact at the reservoir walls. In addition to varying the size and stiffness of 
the tethering potential, we also varied the lattice structure of the tether sites, but settled 
on the simple square and cubic lattices when the pressure-tensor results proved to be 
insensitive to lattice type. 

Systems with Newtonian strips 20 atoms wide were already large enough that no 
systematic deviation from linear stress behavior (with vanishing normal stress difference. 
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a XX = <^yy — 0), resulted. Accordingly we use a narrower system width here, 10, for an 
analysis of boundary driven flows. The system size is iV = 4 x (10 x 10) = 400 particles. 

Time-averaged boundary-driven velocity and normal-stress profiles are shown in the 
figures. The time-reversible frictional forces within the two steadily-moving "boundary" 
chambers were chosen to maintain a kinetic temperature in those chambers (relative to 
the tether velocity) of 0.70. This approximately reproduces the conditions of the Sllod 
and Doll's states in Table II. Fig. 7 shows the average velocity {v^^y)), computed using 
the one-dimensional weight function of Sec. VII. For each of 400 equally-spaced gridpoints 
the instantaneous values of the laboratory-frame velocity components {xj} were spatially 
averaged at the grid points: 

Vxivd) = ^XiWiD{\yG - yi{t)\)/^wiD{\yG - yi{t)\) 

i i 

These instantaneous gridpoint averages were then themselves averaged over time and are 
plotted in Fig. 7. 

Pressure-tensor components at each particle were calculated in two different ways. 
Using the assumed analytic form, a linear velocity profile, 

< y < 10 > (Wx)analytic = +5 ; 

10 < ?/ < 20 > (fx)analytic = 15 - y ] 

20 < ?/ < 30 ^ (t^z) analytic = "5 ; 

30 < y < 40 > (Wx) analytic = 2/ - 35 , 

the pressure tensors for particles in the vicinity of each grid point were computed as 
averages, using the two-dimensional weight function W2d{i^ < h = 3). These individual 
pressure-tensor components were then averaged, with WiD{r < h = 3) as were the veloc- 
ities of Fig. 7. Pressure-tensor components were also computed and averaged using the 
smooth-particle weights W2d{i^ < 3) to calculate the instantaneous flow velocity {v{rj)) at 
each particle j rather than using the assumed linear profile. The differences are relatively 
small, as can be seen in Fig. 8, where the two approaches are compared. Using the in- 
stantaneous velocity at each particle is analogous to, but smoother than, the "unbiased" 



procedure discussed by Evans and Morriss IGj 



The normal stress differences following either approach are considerably larger than 
the Sllod results (Table II), and have the opposite sign to the Doll's homogeneous results 
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Figure 7 
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Figure 7: Velocity profile for a 400-particle system with the aspect ratio illustrated in Fig. 
4. The tethering potential's force constant, k = 100 provides an efficient coupling between the 
driving chambers at < y < 10 and 20 < y < 30 and the driven Newtonian chambers at 
10 < y < 20 and 30 < y < 40. The locations of the two Newtonian chambers are emphasized 
in the velocity plot. The averaged velocity profile was calculated with the smooth-particle 
weighting function wioir < 3). The measured strain rate in the Newtonian regions is about 
±0.48. 

(Table II). 

P — P 

~ 0.03 (Boundary-Driven) . 

The SUod algorithm predicts a smaller effect (smaller by an order of magnitude) while 
the Doll's algorithm predicts the wrong sign! We conclude that the two homogeneous 
algorithms provide no more than an order of magnitude estimate of the normal stress 
effects and further that these effects can be otherwise measured rehably, but with some 
difficulty. 

Fig. 8 shows stresses for the two portions of the four-chamber system with 13 < y < 17 
and 33 < y < 37. These portions have "typical" bulk fluid averages, without any influence 
from the two driving boundary regions. This is a consequence of the smooth-particle 
weight functions' range, /i = 3. 
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Figure 8 
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Figure 8: Boundary-driven two-dimensional flow using the "soft-disk" potential. Average nor- 
mal stress differences are shown. At the left the stresses are calculated relative to the instan- 
taneous velocity profile. There the spatially-averaged instantaneous velocity and stress at each 
particle are computed with W2d, then averaged to get instantaneous profiles using wid, and 
finally time averaged. At the right the stresses are calculated relative to an assumed linear 
velocity profile. In both cases the stress difference is shown for the two regions 13 < y < 17 and 
33 < y < 37 free of boundary influences and hence typical of bulk fluid. The run length was 
5000. 

VIII. NUMERICAL RESULTS: THREE DIMENSIONS 
A. Periodic Shear 

Three-dimensional boundary-driven simulations require only the addition of z coordi- 
nates, with periodic boundary conditions in the z direction. For comparison purposes, 
we first generated series of isoenergetic SUod and Doll's periodic shears. These results, 
shown in Table III, are for periodic shearing of L x L x L cubes of soft- sphere fluid at unit 

density and energy: 

= 100(1 - rY ; = + + ; N/V = Nm/V ^E^K + ^ = 1 . 

The constant-energy ergostat forces keep the total energy of the N — L x L x L particles 
fixed. The kinetic energy X is a sum in which each particle's contribution is measured 
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Figure 9 
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Figure 9: Time-averaged velocity profile for a 4000-particle three-dimensional system with the 
aspect ratio illustrated in Fig. 4. The tethering potential's force constant, k = 100 provides an 
efficient coupling between the driving chambers at < y < 10 and 20 < y < 30 and the driven 
Newtonian chambers at 10 < y < 20 and 30 < y < 40. The locations of the two Newtonian 
chambers are emphasized in the velocity plot. The averaged velocity profile was calculated with 
the smooth-particle weighting function wiD{r < 3). The measured strain rate in the Newtonian 
regions is about ±0.47. The run length was 1000. 

relative to the local velocity. Choosing the cube center as the coordinate origin, the 
systematic velocity in the x direction is taken to be proportional to y: 

{ p^/2m = \pl+pl+pl]/2m ; (p^/m) = v^^ - ey } . 

The pressure-tensor results, given in Table III, are very insensitive to system size. Notice 
that the average kinetic temperature for all these isoenergetic simulations is approximately 
0.5. The Sllod algorithm gives 

Txx > Tyy > (Sllod) , 
while the DoU's-Tensor algorithm gives instead 

Tyy > Txx > %Z (Doll's) . 

The (correct) "Boundary-Driven" results, described next, show instead the ordering 

Txx > Tzz > Tyy (Boundary — Driven) . 

The boundary-driven results also show qualitative differences from the Sllod and Doll's 
results in the normal-stress differences, {Pxx — Pyy) and {Pxx — Pzz) ■ 
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B. Boundary-Driven Shear 



We next implemented 4000 = 4 x (10 x 10 x 10)-particle boundary-driven shear flows 
similar to the two-dimensional flows of Sec. VllB, but with periodic boundaries in the 
z direction and with a flxed isothermal boundary temperature 0.5, chosen to match the 
homogeneous periodic results. For comparison with these Sllod and Doll's results the 
same nominal strain rate, dvx/dy ~ 0.5 was used. Thus the two thermostated chambers 
move with velocities {vx — ±2.5). 

The time-and-spatially-averaged velocity proflle computed with the weight function 
WiD{\5y\ < 3) is shown in Figure 9. The measured strain rates in the straight-line portions 
of the profile are about ±0.47. The energy dissipation rate for 6000 thermostated degrees 
of freedom in the two moving 1000-particle reservoirs was 6000kT{() — 6000x0.5x0.112 = 
336, giving an estimate for the viscosity: 

V = 7^4xternai/(V^e') = 336/(2000 X 0.5^) = 0.672 , 

within two percent of the periodic result, 0.688 from Table 111. The Newtonian shear 
stresses from the 4000 = 4 x (10 x 10 x 10)-particle simulation are ±0.37 in the bulk 
Newtonian regions, corresponding to 

V = -Pxy/{dvx/dy) = 0.74 . 

The time-averaged spatially-smoothed normal stress differences, 

{Pxx - Pyy)l2 ~ 0.01 ; {Pxx - Pzz)/2 ~ 0.01 , 

are only different with marginal significance, and are shown in Figs. 10 and 11. We 
considered three different system sizes, to ensure that these results are insensitive to 
small geometrical changes. Data for = 4 x 10^ (run length 1000), = 4 x 12'^ (two 
runs of length 500; results from only one of them are shown here as the difference between 
the two was insignificant), and = 4 x 14^ (run length 400) are included in the figures. 
Just as before, the Newtonian regions for which the data are plotted are those free of any 
boundary infiuences in the stress averaging. The corresponding Sllod and Doll's values 
for the stress differences, -1-0.003 and -0.016, respectively, are quite different, just as in 
two dimensions. The disparity shows that neither homogeneous algorithm is even close 
to "correct". The actual difference between Pyy and Pzz is apparently quite small, while 

27 



both the Doll's and the Shod algorithms indicate a relatively large difference of order 
±0.04. The statistical fluctuations in the boundary-driven simulations are not quite so 
large as to mask the ordering of the two normal stress differences, 

P — P > P — P 

Somewhat faster/larger computers could make this conclusion more convincing. 



Table III. Soft-sphere viscosities in three dimensions with periodic boundary condi- 
tions. Space-and-time-averaged pressure tensors for homogeneous Shod and Doll's algo- 



rithms using the soft-sphere pair potential illustrated in Fig. 6, 



lOOd 



The 



energy and density are equal to unity and the strain rate, e = dv^/ dy is 0.50. The pressure- 
tensor components, are given in the order xx, yy, zz, xy with the kinetic, potential, and 
total terms indicated. The boundary conditions are periodic, with a total run time of 
200x10,000 timesteps for = 10 x 10 x 10 = 1000 with a fourth-order Runge-Kutta 
timestep of 0.005. 



pK p9 pr. 



pK p<P pl^ 

yy yy yy 
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zz zz 



dK 

xy 



3* 

xy 



xy 



216S 

216D 

512S 

512D 

lOOOS 

lOOOD 

1728S 

1728D 

2744S 

2744D 



0.506 2 
0.496 2 
0.506 2 
0.496 2 
0.507 2 
0.496 2 
0.507 2 
0.496 2 
0.507 2 
0.496 2 



.012 2.518 
.004 2.499 
.010 2.516 
.001 2.496 
.009 2.516 
.000 2.496 
.009 2.516 
.001 2.497 
.009 2.516 
.000 2.496 



0.497 2 
0.507 2 
0.497 2 
0.507 2 
0.497 2 
0.507 2 
0.497 2 
0.508 2 
0.497 2 
0.508 2 



.014 2.511 
.022 2.529 
.012 2.509 
.021 2.528 
.012 2.508 
.021 2.528 
.012 2.509 
.020 2.528 
.012 2.509 
.020 2.528 



0.493 1. 
0.493 1, 
0.493 1. 
0.493 1. 
0.493 1. 
0.493 1. 
0.493 1. 
0.493 1. 
0.493 1. 
0.493 1. 



991 2.483 
991 2.484 

989 2.483 

990 2.483 
989 2.482 
989 2.482 

988 2.481 

989 2.481 
989 2.481 
989 2.482 



-0.062 
-0.062 
-0.063 
-0.063 
-0.063 
-0.063 
-0.063 
-0.063 
-0.063 
-0.063 



•0.282 
•0.281 
•0.280 
•0.281 
•0.281 
•0.280 
•0.281 
•0.281 
•0.281 
•0.280 



-0.344 
-0.342 
-0.343 
-0.343 
-0.344 
-0.343 
-0.344 
-0.344 
-0.344 
-0.343 



As the system size is increased, with the strainrate fixed, the Newtonian temperature 
increases also, in rough accord with the linear model treatment of Sec. II. That is, the 
central temperature increase is proportional to L^. Figure 12 shows temperature profiles 
for the three system sizes considered here. We consider the kinetic temperature here, 
because of its relative conceptual simplicity and its physical conceptual basis 37|]. In the 
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Figure 10 
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Figure 10: Normal stress differences for boundary-driven three-dimensional flows using the 
"soft-sphere" potential. The normal stress difference {Pxx — Pyy)^/'^^ is shown here for three 
different system sizes, with stress calculated relative to the instantaneous velocity profile. The 
spatially-averaged instantaneous particle values of velocity and stress are computed with w^d- 
Then the particle values are averaged to get instantaneous profiles using wid- The figure shows 
time averages of those instantaneous profiles. Only data from the two distinct regions free of 
boundary averaging influences are shown here. 

boundary-driven shear flows the ordering of the kinetic temperatures is 

T > T > T 

with the difference between T^x and T^^ two or three times larger than that between T^z 
and T„„. 




IX. SUMMARY AND CONCLUSIONS 



We were able to characterize the tensor temperature and the nonlinear stresses for both 
homogeneous and boundary-driven versions of simple shear. Neither Sllod nor Doll's gives 
the correct ordering of the kinetic temperatures { Tjj }. Generally the "Sllod" algorithm 
gives a somewhat "better" approximation to the normal-stress differences, { Pa — Pjj }, 
though Sllod is certainly far from "correct" . Despite its evident failures, there is a fairly 



widespread faith in the Sllod approachj60]. It is clear (as emphasized to us by Jim Lutsko; 
see also ref. [61]) that the algorithms' extra rotational terms in the motion equations, 
—epx for Sllod and —tpy for Doll's, when left to their own devices, would eventually cause 
pI and to diverge for the Sllod algorithm, and Py and Tyy to diverge for Doll's. This 
provides a clear explanation of the qualitative difference between the two algorithms' 
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Figure 11: Normal stress differences for boundary-driven three-dimensional flows using the 
"soft-sphere" potential. The normal stress difference {Pxx — Pzz)y/'^^ is shown here for three 
different system sizes, with stress calculated relative to the instantaneous velocity profile. The 
spatially-averaged instantaneous particle values of velocity and stress are computed with w^d- 
Then the particle values are averaged to get instantaneous profiles using wid- The figure shows 
time averages of those instantaneous profiles. Only data from the two distinct regions free of 
boundary averaging influences are shown here. 

predictions and the corresponding opposite directions for the rotation of the principal 
axis of the stress. 

Lutskoj62| has reviewed the hard-sphere-based Enskog theory for nonhnear stress ("the 
only viable theory") and his finding that Pxx > Pyy in simple shear is quite consistent 



with our results. On the other hand some theoretical models [63l| and some computer 
simulations jl6 1 find Pxx < Pyy, even for relatively simple fluids, so it is clear that more 
investigations are required. Evidently the temperature tensor and the nonlinear stresses 
are not given accurately by the Shod algorithm. The more realistic boundary- driven flows 
need to be used whenever confidence in the results is required. 

Boundary-driven flows are actually extremely complex, even for this simplest possible 
model of shear. The flows we can study are dominated by fluctuations which can be tamed 
by averaging, in one, two, or three dimensions, but the time-averaged flows describe the 
time-dependent situation no better than they would for a physical waterfall or a turbulent 
stream. 

It is fortunate that a hydrodynamic description of flows is feasible on a very small scale 
(just a few particle diameters) as was apparent from the earliest shockwave simulations, 
which showed shockwidths of only a few particle diameters. It still remains a puzzle 
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Figure 12: Tensor temperatures for boundary-driven three-dimensional flows using the "soft- 
sphere" potential. The time-averaged kinetic temperatures, relative to the instantaneous velocity 
profile, computed with w^d and wid, are shown here for three system sizes, all with boundary 
temperatures and strain rates equal to 0.5. The systems contain four chambers with 10^ (run 
length 1000), 12 (run length 500), and 14^ ( run length 400) particles per chamber in the three 
cases shown. Although data from the two distinct regions free of boundary averaging influences 
are shown here the differences between Txx, Tyy, and T^z, are of order ±0.02 and are too small 
to see on the scale of the figure. 



that Shockwaves indicate an enhanced nonhnear viscosity while the homogeneous shear 
algorithms considered here predict a reduction rather than an enhancement 45]. 
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